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Q3l)G - A COMPUTER PROGRAM FOR STRA \ N-ENKRGY-RELKASE RATES FOR 
OKI. AM I NAT I ON GROWTH I N COMPOS ITK LAMINATES 


FORTRAN V APPROX I MATELY 4000 SOURCK STATKMKNTS 
9 TRACK 1600 BPt CI)C NOS INTERNAL FORMAT MAGNETIC TAPE 

ABSTRACT 

The Q3l)G is a computer program developed to perform a quasi-three- 
dimensional stress analysis for composite laminates which may contain 
del ami nat Ions • The laminates may be subjected to mechanical, thermal and 
hygroscopic loads. The program uses the finite element method and models the 
laminates with eight-noded parabolic isoparametric elements. The program 
computes the strain-energy-release components and the total strain-energy 
release in all throe modes for delamination growth. A rectangular mesh and 
data file generator, DATGEN, is included. The DATGEN program can be executed 
interactively and Is user friendly. The doeumentat ion includes sections 
dealing, with the Q3D analysis theory, derivation of element stiffness matrices 
and consistent load vectors for the parabolic element. Several sample 
problems with the input for Q3DC and output from the program are included. 

The capabilities of the-v DATGEN program are illustrated with examples of 
Interactive sessions. A microfiche containing all the examples presented in 
this report is included with the documentation. The Q3DG and DATGEN program 
have been Implemented' on CYBER 170 class computers. Q3DG and DATGEN wore 
developed at the Langley Research- Center during the early eighties and docu- 
ment ed in I 984 - 1 983 . 
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INTRODUCTION 


Because of their high s t rength-to-we ight ratios, composite materials are 
being extensively used in aerospace, mechanical, and automotive structures. 

As these materials are being considered for heavily loaded primary aircraft 
structures, increased attention is being focused on understanding and charac- 
terizing composite do 1 am i na 1 1 on * Del ami nat ions may result from low-velocity 
impact, eccentricities in load paths that induce out-of-plane loads, or 
discontinuities in the struct urc* which create out-of -plane stresses. 

Fracture mechanics concepts applied to delaminated structural components 
wore found to be useful in understanding and characterizing delaminations. 
Strain-energy-release rates associated with delamination growth were shown to 
be useful in predicting the onset of edge del ami nat ions in composites [l, 2]. 
While the total strain-energy-release rates can he calculated using the 
classical laminated plate theory (CLT), the separation of the total strain- 
one rgy- re 1 ease rates into the three components associated with the three modes 
of fracture cannot be accomplished using the CLT approach. To calculate the 
individual components, one has to resort to numerical methods such as finite 
element analyses. 

The purpose of this report, therefore, is to document the computer 
program, Q3DC, which uses a finite element analysis to calculate the strain- 
energy- re lease rates for long, rectangular composite laminates containing 
de 1 ami nat i ons . The laminates may be subjected to mechanical, thermal and 
hygroscopic loadings. The Q3DG program uses the quas i-three-dlmens ional (Q3D) 
analysis in conjunction with linear fracture mechanics concepts and finite 
element methods. The program calculates , Gjp and Grjr, the strain-energy- 
re lease rate components due to tin? throe fracture modes, and the total strai n - 
energy-release rale (I for de lamination growth in composite laminates. The 
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program can also be used as a Q3D stress analysLs tool for laminates with or 
without de laminations . 

First, the theoretical aspects of the Q3D analysis arc presented. Next, 
the finite element analysis and the st ra I n-energy-rel ease rates are discussed. 
Then, the program specifications and organization an* presented. Finally, the 
input data for and output of the program are described. Appendix A presents 
the derivation of the stiffness matrix and thermal vectors for an eight -noded 
isoparametric element. Appendix B lists the names and functions of subrou- 
tines and major program variables and the common blocks and their elements. 

In Appendix C, three example problems with their output are presented. 

Appendix D illustrates some conditions which cause? the program to terminate 
and discusses some debug strategies. Next, Appendix H describes DATGKN, an 
automatic rectangular mesh generator, and shows several example problems using 
the mesh generator. Appendix F presents the COSMIC documentation require- 
ments. Finally, Appendix G describes the floppy disks for the IBM PC and 
compatibles. These disks can be easily reproduced and used to upload the 
source and data files to the mainframe computers. 

DESCRIPTION OK THE DF.LAM IN AT [ON PROBLEM 
Consider a long, rectangular laminate will) an edge del am (nation as shown 
in fig. 1. Away from the ends where the loads are applied, the displacements 
at any x = constant plane are assumed to be given by (refs. 4-8) 

u ( x , y, z) = e o x + U(y, z) 

v ( x , y , z) = V(y , z) ( ! ) 

w(x, 'y , z) = W(y , z) 

The term is the uniform axial strain and IJ, V, and W are functions 

of y and z alone. Equation (l) describes a "quas i-three-d i mens i onal " 

(Q3D) problem. The modifier "quasi." is used because there are displacements 
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in l hn*«* d I reel Ion?:, bul l li«» gradients of II, V, and W with respect to the 
x-coord ( nat e are zero. 

Because of the symmetries inherent in the layups of the fiber reinforced 
composites, the displacement functions U, V, and W satisfy the following 
requirements (4-7]. 


U(y, z) = -U(-y, -z) 

V(y , z) » -V(-y, z) (2) 

W(y, z) = -W( y, -z) 

For mechanical loading, the uniform axial strain c o in equation (1) is 
specified. For thermal and hygroscopic loads the laminate Is subjected to 
temperature and moisture changes, AT and AM, respectively. For these cases 
the magnitude 1 of e 0 in equation ( 1 ) is unknown and is determined as a part 
of the solution (see ref. 9, 10 ). 

Tn the Q3!) analyses, each ply is idealized as a homogenous, elastic 
orthotropic material. The material properties of each ply are defined by 
Young’s moduli i (Ej|, ^2’ £ 33 ), shear moduli! (G^ 2 > ^13» ^ 23 ), Poisson's 
ratios (v| 2 > v \]> ^ 3 ), and ply fiber angle (0). In addition, the expan- 
sion coef fi dents for thermal and hygroscopic loadings, U| , a 2 , and 03 

and [ij , 02 * and ^ 3 , respectively, are needed. The subscripts, l, 2, and 3 

o n the various quantities listed above correspond to the longitudinal, trans- 
verse, and thickness directions, respectively, of a zero-degree ply* 

Finite Element Analysis 

Because of the symmetries In the problem, only one quadrant of the 
x - constant plane is idealized. This quadrant (0 < y < b and 0 < z ^ t) is 
shown as the 1 shaded region in fig- 1(b). The analysis region is modeled by 
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eight-noded , isoparametric elements fill. At the del ami nat Ion tip, no singu- 
larity elements are used. Instead, the finite' element Ideal i za ( I on Is usually 
made much finer near the delamination tip. A typical finite clement ideal- 
ization involving 367 nodes and 102 elements is shown in fig. 2. Due to 
symmetry the displacement functions, U and V, are prescribed zero at all 
nodes on the y = 0 line, and the d i spl acoment function W is proscribed zero 
on the z ~ 0 line. 

The derivation of stiffness matrices and the mechanical, thermal and 
hygroscopic load vectors needed for the finite element analysis is described 
in Appendix A. 


St ra i n-Kne rgy-Re lease* Rates 

Mode 1 (opening mode?), mode* TI (sliding mode), mode 111 (tearing mode) 
components and the total strain-energy-release rates for del ami nation growth 
are calculated by a virtual crack extension technique similar to that reported 
in ref. 12. This technique uses the nodal forceps transmitted at and near the 
delamination tip and the relative displacements just behind the del ami nation 
tip to determine the strain-energy-release rates. The nodes at which the 
forces and relative displacements are evaluated are shown In fig. 3. The 
s tra Ln-energy-re lease rates are computed as shown hi' low: 


C = fp (W, - W. ) + I- (W - W )|/(2A) 
l z . 1 k z . n in 1 

1 1 


G tt = K (V, - V ) + F (V - V ) 1/(2 A) 
it L y.l k y.n m 1 

1 .1 

G rrr = I F < V , ' V + F < u ~ ^ > |/<2A> 

III 1 x t 1 k x.n m * 

* i 


G 


G i + G n 


+ G 


rn 


O) 


A 



rho forces b [• and I* , are the To roes in x-, y-, and 

x , y , x . 

f i i 

z-d i rec t i ons , respect i vol y , acting at node 1. The forces at node i, shown 
tn fig. 3, are computed from the elements A and B while the forces at node 
) are computed from element A alone. This procedure for computing the 
strain-energy-release rates requires that the finite element mesh he symmetric 
about the crack tip. Square elements with a side of length A - h/4 were 
found to be the optimum size and, hence, an arrangement of the type shown in 
fig. 3 is recommended* 

The program first computes the forces and displacements individually 
for a mechanical strain of one micro in/ in (m/m), a temperature change of 
+ i°F (°K), and a moisture change of one percent. Then, the program computes 
the appropriate forces and displacements for the given input values of mechan- 
ical strain and temperature and moisture changes and uses the forces and dis- 
placements In equation (3) to calculate the s t ra i n-energy-rel ease rates. The 
mode I, r 1 , Ilf components, and the total strain-energy-release rates are 
computed for the seven possible loading combinations listed below: 

1. Mechanical loading acting alone. 

2. Thermal loading acting alone. 

3. Hygroscopic loading acting alone. 

4. Mechanical and thermal loading acting together. 

r >. Mechanical and hygroscopic loading acting together. 

b . Thermal and hygroscopic loading acting together. 

/. Mechanical, thermal, and hygroscopic loading acting together. 

PROGRAM SPKC l F I CAT IONS 

The program Q3l)G is writttui in FORTRAN V and has a core requirement of 
about 343,000 octal words. The program was compiled and successfully executed 
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on CDC CYBER 170 class computers. The current maximum allowable values of 
major arrays are given in Appendix B. 

The program uses n scratch unit. Tape 10, for storage of the element 
stiffness matrices and consistent load vectors when they are first computed. 
This information is later used to check the element and global equilibrium 
conditions. A scratch unit is used since the recomputation of element stiff- 
nesses and load vectors Is more time-consuming I ban reading from L he scratch 
unit. 

-r i 

The program uses an In-core equation solver. Hence the maximum size of 
the finite element model is limited by the memory of the computer being used. 
Most of the core is required to hold the global stiffness matrix, B1GK. The 
array BIGK is dimensioned (7b, 1110), which permits l I 10 degrees of freedom 
(370 nodes) and a bandwidth of 75. 

The Q3DG program can be easily modified using a text editor to change any 
of the problem dimensions. The required changes are given below: 

1) Change the string " 75,1110 11 to " xx,yyy " everywhere, where xx and 
yyy are the new number of rows and columns, respectively, of the array BIGK. 

2) Change the string (83250) to "zz^” everywhere, where /./z is the 
product of xx and yyy. 

3) Change the string "(1 110,3)" to "(yyy, 3)" everywhere, where yyy is 
the number of columns in the array 8IGK. 

4) Change the string "X(400, 2), N0I)(200, 8)" to "X(aaa, 2), N0i)(bbh, 8)" 
everywhere, where aaa is the maximum number of nodes and bbb is the maximum 
number of elements. 

5) Change the string "(400)" to "(aaa)" everywhere, where ana is, again, 
the maximum number of nodes. 
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6) Change the string "(200)" to "(bbb)" everywhere, where bbb, is, 
again, the maximum number of elements. 

7) Change the DATA l.X, . . . card which is the list line in the source 
code. The Integers, LX, LNOD, LUST, LOU), are computed as follows: 

LX = aaa * 2 

LNOD = bbb * 8 

LUST = 200 

LOLD - aaa 

8) Change the dimensions of the arrays NDUM, NDUI), and XDUM in 
subprogram ADJUST to the following: 

NUUM( bbb, 8) , NDUD(bbb) 

XI)UM(aaa,2) 

and change the equivalence string ", (BIGK (1601), NDUD (1))" to ", (BIGK 
(bbb*8+l), NDUD (l))". 

9) Change the string "0(1 110)", the working storage for the band- 
solver, in subprogram ASEKBL, to "D(yyy)". 

10) Change four lines In subprogram FORCES as follows. 

DIMENSION FORCE (aaa, 3,3), SP(aaa, 6,10), ND(aaa,10) 

EQUIVALENCE (BIGK(I), FORCE (1,1,1)), (BIGK (9*aaa+l), SP(l,l,l)), 

(BIGK (aaa*9Faaa*60+l), ND( 1,1)) 

DIMENSION Si'A (aaa, 6), NDA (aaa) 

EQUIVALENCE (BIGK (aaa*9 f aaa*360 +• aaa* 10 + 1), 

SPA (I)), (BIGK (ana*9 F aan*360 + aaa* 10 F aaa*6 F 1), NDA (1)). 

The program is set up for a 3 * 3 Gaussian integration of the element 
stiffners matrices. The order of integration can be changed by changing the 
f o I I ow i ng, : 
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L. NGAUSS in the BLOCK DATA subprogram (currently it rends 
DATA NGAUSS J 3 1 

change 3 to n, where n is the order of integration desired (2 < n < 8). 

2. If the order of integration n is chosen to be greater than 3, 
change the dimensions of array PDER everywhere in the program to (n*n*2,8). 
Current dimension of array PDKR are (18,8). Also note that I lie array PDER is 
the COMMON block ODER. 

PROGRAM ORGANIZATION 

In this section the flow of Q3DG is described. Brief descriptions of the 
subprograms and major program variables are presented in Appendix B. 

Fig. 4 presents the flow chart of the program. Most of the input Is read 
in the main program Q3DG, which calls 5 main subprograms. Subprogram MODULUS 
generates the modulus matrices and their thermal and hygroscopic strain 
vectors for all the materials. Subprogram CIJ) calculates the coordinates of 
the midside nodes (if these are unspecified), roads boundary conditions and 
creates the boundary condition array, LIST. Subprogram ADJUST reorders the 
coordinate array and the nodal connectivity array according to the renumbering 
scheme Input by the user. Subprogram ASKMBL calls the relevant routines to 
evaluate the stiffness matrices and load vectors and assembles the global 
stiffness matrix. Subprogram FORCES checks equilibrium at the element and 
global levels, calculates the nodal stresses and calls the G-cal cul at Ion 
routines. The functions of the subprograms called by ASEMBL and FORCES are 
explained in Appendix B. 


8 



INPUT DATA 


The required input data is described in this section. The data can be 
created on a file, LFN, equated to TAPES. An alternate method of creating 
the input file is to use the data generation program DATGEN described in 
Appendix E. 


CARD 

NUMBER 

COLS 

FORMAT 

VARIABLE 

DESCRIPTION 

SET 

OF CARDS 





1 

1 

1-80 

20A4 

TITLE 

TITLE of problem 

2 

1 

1- 8 

A8 

POUT 

Output option. SHORT — 
for short output. XLONG 
for long output. 


1 

1- 5 

r r > 

NPOIN 

Total number of nodes 



ft- 10 

i r > 

NELEM 

Number of Elements 



11-15 

15 

NFRKE 

Number of degrees of 
Freedom per node = 3 



16-20 

15 

NNODE 

Number of Nodes per 
element = 8 

4 

* 

1-10 

F10.5 

Y1 

Coordi nate 



1 1-15 

15 

JC0RD( 1 ) 

Node numbers with this 



16-20 

15 

JC0RD( 2) 

coordinate in columns 



• 

♦ 

• 

11 through 60. 



• 

56-60 

• 

15 

• 

JCORD( 10) 



2 sets of cards of this 
format, with each set 
terminated by a zero (or 
blanks in cols. 11-15). 

First set is y-coordi- 
nates; second set is 
z-coordinates. 

* l nput until nil y coord i nates are specified. End y-coord i nates with a blank 
card. 

input, until all z~coord i nates are specified. End z-coord i nates with a blank 
card. 


5 

NELEM 

I- 5 

15 

l 

Element number 



6- 10 

15 

N0D(I , l) 

Element connectivity 



11-15 

15 

NOI)( 1,1) 

starting at any corner 


node in the counter 
clockwise direction. 


. • 

[5 NOI)( I , NNODE ) 

One card for each element 
hence, NELEM cards. 
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CARD 

NUMBER 

COLS 

FORMAT 

VARIABLE 

DKSCR I t*T 1 ON 

j>ET 

OF CARDS 





6 

i 

l- 5 

15 

NMAT 

Number of materials 

7 

4 *NMAT 

1-20 

E20.7 

EM( 1,1) 

Specification matrix for 



21-40 

• 

E20.7 

• 

EM( 1,2) 

• 

materials. 




• 

E20.7 

• 

F.M(NMAT, 16) 

This card occurs A times 
for each material (NMAT 
materials), containing A 
fields per card: total of 
16 el eme n t s pe r material. 

8 

t 

1- 5 

15 

NF 

Ply-Material correlation. 
One card per ply. 

Number of first element 
in ply 



6-10 

15 

NI, 

Number of I a s 1 e 1 e me n t 



11-15 

15 

N! 

increment 



16-20 

15 

MR 

Material number 


e.g. 

1 

7 2 

K 

Defines elements 1, 3, 5, 
and 7 to belong to the 
material . 

tRepeat until 

all elements are defined. Terminate this set of cards 



with a 

zero (or 

blanks) In cols. 1-5. 
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1- 5 

15 

NPLY Number of plies In 
the laminate 




6-10 

15 

iRHS Number of loading 
conditions (see card 10). 




1 1-20 

F10.4 

WIDTH Width of lamina to 

10 

i 

1-20 

E20.7 

SMKOH 

Meehan 1 oa 1 St ra i n 



21-40 

E20.7 

DKLT 

Temperature change 



41-60 

E20.7 

DKLH 

Hygroscopic (moisture) change 

11 

i 

1- 5 

15 

LEM( l) 

Elements used In the G- 



6-10 

15 

LEM(2) 

calculation (elements A fit B 
shown in Fig. 3.) 



e-g. 

44 

45 

For the model shown in fig. 5 

12 

i 

1- 5 

15 

NGF( 1) 

Nodes at which forces are 



6-10 

15 

NGF(2 ) 

evaluated and used In the G- 
calculation (nodes i, j in 
Fig. 3 and in that order.) 



e.g. 

205 

196 

For the model shown in fig. 5 
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CARD 

NUMBER 

COLS 

FORMAT 

VAR I. ABLE 

DESCRL PT ION 

SET 

OK CARDS 





13 

i 

1- 5 

15 

NGD(l) 

Nodes at which relative dis- 



6-10 

15 

NGD(2) 

placement are evaluated and 



11-15 

15 

NGD(3) 

used in the G-calculat ion 



16-20 

15 

NGD(4) 

(nodes k, 1, m, n in fig. 3 and 
in that order). 



e.g. 

226 227 

216 217 

For the model shown in fig. 5 

14 

t 

1- 5 

15 

NN 

Node Number 



6-10 

15 

NX 

Restraint code in the X- 
di reel ion 



11-13 

15 

NY 

Restraint code in the Y- 
d irect Lon 



16-20 

15 

NX 

Restraint code in the Z- 
direction. One card for each 
node with boundary conditions. 
Terminate this set of cards 
with a zero (or blanks) in 
cols. 1-5. 


Rest ra I n t 

Code: 

0 Is f ree 2 

l is fixed 


tRepeat until all boundary conditions are defined. 

IS I i- 5 15 IRFNUM Renumbering option: 0 if 

no renumbering desired; 1 
if renumbering scheme 
fol lows. 

lb * 1-80 1615 JNKW(NPOIN) 

Renumbering scheme: 
needed only if IRENUM not 
equal to zero 


* 


NP0IN/I6 (rounded up to the next integer when necessary) 



Sp ecia 1 C a se s_ 

a ) De l a mi na t ion a t m i dp_l a ne ( z = 0 ) 

When the delamination Is at the midplant' and the laminate Is symmetric, 
the input is slightly different* 

Card sets 1-10 and card sets 14-16 remain unchanged* The O-cal culat i on 
card sets 11-13 need to be input slightly differently. To explain this, con- 
sider the modeling near the del ami nation as shown in fig. r >(h). Card sets 
11-13 will be as follows: 


Card 

set 

ii. 

44 

45 

Element used in the G-calcul at ton . 

Card 

set 

12. 

205 

196 

Nodes at which the forces are 






evaluated and used in the G- 






cal cul at ions 

Card 

set 

13. 

-227 227 

-217 217 

Nodr's at which relative displace 


meats are evaluated and list'd in 
the G-cal culat I ons . 


Note the negative signs on the integers NGD(l) and NG1) (3). The program rec- 
ognizes that deformation is symmetric about the mi dpi ant' and calculates the 
appropriate relative displacements at these nodes. A 1st) nol r that for this 
case, only mode I st ra 1 n-ene rgy-re 1 ease rale exists, intuit » II ami moth* 111 art* 
identically zero. 

b) Quasi-thr ee-dim ensi onal analysis w i t h o r wi t hou t 

The program can be used to perform edge-stress analyses very easily for 
laminates without delaminat ions and stress analysis alone (without 
G-calculat ions) for laminates with del ami nat ions . To exercise this option, 
simply input zeros in card sets 11, 12, and 13 as shown below: 

Card set 11. 0 0 

H " 12. 0 0 

" " 13. 0 0 0 0 

Note that all three card sets 11, 12 and 13 must be input. 
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OUTPUT 


The output of the Q3DG Program Ls presented In this section. Two ver- 
sions of output, long and short, are possible with this program. If the user 
exercises the long output option (XLONG on the second card in columns 1-5), 
then the following output will be obtained. (Only the general description of 
major output sections Ls presented.) 

*T L t le 

^Output option 
^Program specifications 

*Y-coord I nates and corresponding node numbers 
*Z-coord Lnates and corresponding node numbers 

*Materlal Specifications: Modulus matrix, thermal and hygroscopic strain 

vectors for all materials used In the model. 

*Ply-material correlation 

*C-calculat ion variables: Elements and nodes involved in the G- 

calcul at ions . 

^Boundary condl t Ions 
*Ronumbe r I ng opt I oils 

Element stiffness matrix (24 x 24) and consistent loads for element 
number l . 

Sum of F x , F , and for all nodes before and after boundary 

condl t Ions. 

^Degrees of freedom and bandwidth: Maximum and current values. 

Displacements (U, V, and W) at all nodes for each of the loading 
condl t ions. 

Nodal forces at all nodes for the elements Involved In the G-calcul ation. 
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*Sum of forces la the model due to each of the loading conditions. 

Strain e x required to satisfy E F x = 0 for thermal and hygro- 
scopic loadings. 

Nodal forces at each of the nodes for each of the loading conditions. (if 


all of the nodal forces, F , F 

X, y, 


and K , at node I are less 
7 -i 


than 1.0E-6, then the forces are not listed). 

Average nodal stresses at each of the nodes In the model due to combined 
loading. 

Average nodal stresses at each of the nodes In each material due to 
combined loading. 

^Equilibrium checks for each of the loading conditions. (EF x , EF^, EF^, 


t 


I 


i 


1 = 1, Node? in the model). 

t*St rain-energy-release rate calculations - 

Forces and displacements involved in the O-calculations for each of 
the loading conditions# 

Gf, Gxi> and for each of the seven possible loading 

combinations. 

t*Summary of the strain-energy-release rates for the seven possible 
loading combinations . 

Element equilibrium. If all elements satLsfy equilibrium, l.e., if 

EF X = EFy = EF Z < 1.0E-6 for each element, then the program prints 
M A1 i elements satisfy equilibrium. The solution may be correct.” 

If some elements do not satLsfy equilibrium the program prints a 
warning message that n numbers of elements do not satisfy 
equilibrium and lists the n element numbers. 
tThese items do not appear in the output when only Q3D stress analysis is 
performed (i.e., when card sets 11, 12, and 13 contain zeros). 
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Only the asterick (*) items appear In the output when the user exercises 


the short output option (SHORT on the second card in columns 1-5 )• 
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APPENDIX A: DERIVATION OF STIFFNESS MATRTX AND CONS [STENT LOAD VECTORS 


This appendix presents the derivation of the stiffness matrix and the 
thermal load vector for the eight-noded element used In the Q3D analysis. The 
derivation of the hygroscopic vector is very similar to that of the thermal 
vector and hence is omitted. Additional details on these derivations can be 
found in refs. 5, 9, and 10. 

As previously pointed out, the Q3D analysis determines the displacements, 
U, V, and W (eq. (1)), at each node In the finite element model. There- 
fore, the displacements at any interior point In the element, 

{u} = (U V W} (A-l) 


are interpolated by the standard shape functions [N ] (Ref. 9) as 

{u} - [N] {6} (A-2) 

where {6} is the vector of element nodal displacements (dimensioned (24,1) 
for the eight-noded element). The strains in the element are obtained by 
dif ferent Lating eq. (A-2) to obtain the following equation: 


U> = IB] {6} 


(A-3) 


where the column vector of the strains (e) is 
{e> = {e x e y e z e xy e yz e zx } T 
Equation (A-3) is partitioned as shown below: 



The strain-displacement relationship Is partitioned fn this 
each element has a prescribed strain of e x = e Q (see eq . 


(A-4) 


( A-5) 

manner because 
( 1 ) ) . The 
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remaining strains (e R ) are the unknown strains in the element and are thus 
related to nodal displacements {6} by the Following equation: 

{e R } - [B r ] {6} (A-6) 

The vector {e R } is of size (5,1) and the matrix [B R ] is dimensioned (5,24) 
for the eight-noded element. 

The stress-strain matrix of the orthotropic material of the element also 
can be written in partitioned form as 



The matrices [D xx ] , [D xR ], an< * [Drr] are of dimensions (1,1), (1,5) and 
(5,5), respectively, and the column vectors {e^}, (e R }» an( * are 

of dimensions (1,1), (5,1), (1,1), and (5,1), respectively. 

The strain energy of the element, then, is written as 

U' - ] {o} T {eJ dvol (A-9) 
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Substituting eqs. (A-8) into (A-9), using the partitioned form of (e), and 
rearranging terms produces the following equation.: 



T 

2e D T 

O XX X 


T T 

2t I) „ T d + T 
o xR R x 


I) T 

XX X 


* 2T X d xR t r + t r d rr V d ' , °' 


+ -i I K °RR C R + 2e R °xR “o 


2 \ '>RR T R 


- 2£ R D xR T x) ^ 


( A- 1 0 ) 


The first integral term In eq. (10) is not a function of the unknown strains 
{e^}, and hence can be disregarded from further consideration. 

The unknown strains are related to the nodal displacements by eq . (A-6). 
Therefore, the strain energy can be written in terms of the nodal displace- 
ments (<$} by substituting eq . (A-6) into eq . (A-10). in the absence of any 
other external loading the strain-energy U* is equal to the total potential 
energy JTp. Minimizing tt^ with respect to the nodal displacements {<5j 
yields the following equation: 

[K]{6| + (f) e + lr > at = 0 ( A- 1 1 ) 

o 
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whore 


ixi / I b r 1 I'Vr 1 1 b r 1 dvo1 

{n e = / [bJ][|^ r ]U 0 } dvol (A- 12) 

o 

tf >AT - -/([»r][»Ir]( T x I * [“r]1 I> RrH T R ) ) -‘"O' 

The matrices [Br], [Orr] , and [D x r] are dimensions (5,24), (5,5) and 

(1,5), respectively; the column vectors {e 0 }, (T x ), and (Tr) are of 

dimensions (1,1), (1,1), and (5,1), respectively; the stiffness matrix [K] and 

the consistent load vectors, {f } £ , and {f) AT are of dimensions (24,24), 

o 

(24,1) and (24,1), respectively. 
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APPENDIX B: SUBROUTINES, MAJOR PROGRAM VARIABLES, AN!) COMMON BLOCKS 


This appendix first presents the names and functions of the subroutines 
and major program variables and, next, gives the common blocks with their 
elements and the subroutines whLch use the common blocks. 


B . I Subrout 1 nes 


Name 


Func t Ion 


1. ADJUST Reorders the nodal coordinate array, the element 

connectivity array and the boundary condition array 
according to the renumbering scheme provided by the user. 

If no renumbering is given, no reordering of these arrays 
is performed. 

2. ASEMBLE Processes all elements, obtains element stiffness matrices, 

load vectors and assembles the global stiffness matrix. 

Also calls the band-solver. 


3. ASTAR 


Obtains the transformation relationship between the 
generalized coordinates and the nodal coordinates. This 
subprogram is called only once, when the first element Ls 
processed. 


4 . BLOCK DATA 


Contains the Gaussian coordinates and weights up to an 
8-point integration rule. Also contains the parent 
coordinates (£,n) of the 8-noded element. 


3. BMAT 


Evaluates the stiffness matrix and the load vectors 
eq. (A-12)) at the integration points. 


(se* 


6. BOUND 


Prescribes the boundary conditions, 


7. CKBNDW 


8. CKDEGF 


9. CLD 


Calculates the current bandwidth from the element nodal 
connectivity and checks if the current bandwidth is greater 
than the maximum bandwidth in the program. 

Calculates the degrees of freedom (dof) in the model and 
checks If the current degrees of freedom are greater than 
the maximum in the program. 

Calculates the coordinates of the midside nodes and reads 
in and forms the boundary condition array. 


10. DERIVE Obtains the derivatives of the shape functions at any point 

(£,n) In the element. 
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1 I . FORCES 
12. GEVAL 

n. int 

14. INTSX 

l 5. MATLNV 

16. MATMUL 

17. MODULUS 

18. PARENT 

19. PICKD 


Calculates the element forces and checks the element 
equilibrium. Calculates the model stresses and forces and 
checks global equilibrium. Calculates the e x strains 
needed to satisfy the condition = 0 for the thermal and 

hygroscopic loadings. 

Computes the strain-energy-release rates for delaraination 
growth for all seven possible loading conditions. 

Obtains the stiffness matrix and load vectors of an eight- 
noded element by Gaussian quadrature. 


Obtains the total force, F x , in the x-direction in an ele- 
ment due to mechanical strain of e x - 10 , a temperature 

change of AT = 1°F, and a moisture change of AH = 1%. F x 
is obtained by Integration of o x stresses on the element. 


Obtains the inverse of a square matrix. . 

Obtains the product of two matrices. 

Computes the modulus matrix, and thermal and hygroscopic 
strains of each material. 


3N. 

Calculates the parent derivatives 
each of the (NGAUSS)^ Gaussian points. 


3N. 

i 


3n 


i - 1,8 at 
Also calculates the 
x 2 Gaussian-nodal stress transformation matrix. 


Calculates the displacement differences (see eq. (3)) 
needed i n the G-cal culations . 


20. PICKF 
2 I . PROCESS 

22. Q3DG 

23. SBANDPl) 

24. SHAPE 

25. SMALLK 


Calculates the forces (see eq.(3)) needed in the 
G-cal cul at ions. 

Summarizes the strain-energy-release rate results for all 
the 7 possi bl e load ing conditions . 

Main program of this analysis. Reads most of the input. 

Solves a linear system of equations Ax = B. A is a banded 
symmetric positive definite coefficient matrix and B is a 
matrix of right hand side vectors. 

Shape functions of the eight-noded element. 

Calculates the element stiffness matrix. 


26. STRESS 

27. TRANS 


Calculates the nodal stresses for an eight-noded element 
from the stresses at the 2x2 Gaussian points. 

Obtains the transpose of a matrix. 


28. ZEROLN Zeros out an integer array. 

29. ZEROLV Zeros out a real variable array. 
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VARIABLE NAME 
AMOD (6,6) 

BIGK (75, 1110) 

CORD (8,8) 

DELH 

DELT 

DIS (1110,3) 

EE (6,6,10) 

ELI) IS (24,3) 
EM (10,16) 

EMO (6,1) 

EPSM 

EPST 


Major Pro gra m Variables 


COMMON 


DESCRIPTION 


s 


MOD 
AST IF 

GAUSS 

CONFG 

CONFG 

DISP 

MOD 


MATER 


MOD 


Modulus matrix of the material being 
processed 

Assembled global stiffness matrix. First 
subscript corresponds to the bandwidth; 
second subscript corresponds to the degrees 
of freedom In the model. 

Gaussian coordinates up to 8-point Gaussian 
Integra 1 1 on . 

Change In moisture content, AH, 

(percentage) of the laminates for 
hygroscopic loading. 

Change in temperature, AT, of the laminate 
for thermal loading. 

Global displacements. First subscript 
corresponds to the degree of freedom; 
second subscript corresponds to the three 
loading conditions. 

Modulus matrices for each material. The 
third subscript corresponds to the material 
number. 

Nodal displacements for l lie three loading 
conditions of tin* element being processed. 

Material property specifications of each 
material. The first subscript corresponds 
to the material number. 

Hygroscopic straLns due to AH - 1% of the 
material being processed. 

Strain e x needed to satisfy the condition 
F x = 0, due to AH = 1%. 

Strain e x needed to satisfy the 
condition F x = 0, due to AT - 1 ° . 
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VARIABLE NAME 
ET (6,1) 

RTM (6,10) 

ETH (6,10) 

FEMT 
FKTT 
EE XT 

FOR (24,3) 

FORCE (400,3,3) 

OCR (8,2) 

COLS (2,3,3) 

GFOR (2,3,3) 


COMMON 


DESCRIPTION 


MOD 

MOD 


MOD 


GENRE 


GCAL 


GCAL 


Thermal strains due to AT = 1 ° for the 
material being processed* 

Hygroscopic strains due to AH * 1° for all 
materials. The second subscript corre- 
sponds to the material number. 

Thermal strains due to AT * 1% for all 
materials. The second subscript corre- 
sponds to the material number. 

Net force in the x-direction of the 
laminate due to AH = 1% and e x = 0. 

Net force in the x-direction of the 
laminate due to AT = 1° and = 0. 

Net force in the x-direction of the 
laminate due to a mechanical strain of 
e x - 10” 6 alone. 

Nodal forces of the element being pro- 
cessed. The second subscript indicates one 
of the three loading conditions. 

Nodal forces at each node of the model. 

The first subscript indicates the node 
number; the second subscript indicates the 
x-, or z-direction; the third subscript 

Indicates one of the loading conditions. 
Array equivalenced with BIGK. 

Parent coordinates (£,n) of the 8-noded 
isoparametric element. Note that 0 < F, 
n < l . 

Relative displacements used in the 
G-cal culat Ion. The first subscript 
Indicates one of the sets of nodes just 
behind the delamination tip; the second 
subscript indicates the x-, y-, or 
z-direction; the third subscript indicates 
one of the three loading conditions. 

Forces used in the G-cal culat ion . The 
first subscript indicates one of the set of 
nodes at and ahead of the delami nation tip; 
the second subscript indicates the x-, y-, 
or z-direction; the third subscript 
Indicates one of the three loading 
conditions. 
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VARIABLE NAME 

COMMON 

DESCRIPTION 

IMAT (200) 

MATER 

Array containing material code (number) of 
all the elements. 

IMT 

4 MOD 

Material code of the element being 
processed . 

f RENUM 


Renumbering option. When rRKNUM=0, no 
renumbering option Is exercised. User must 
provide the renumbering scheme. 

IRHS 

CONFG 

Number of loading conditions. Maximum 
number' of loading conditions that can be 
specified Is 3: mechanical, thermal and 

hygroscopic loadings. 

JNEW (400) 

RENUM 

Array which relates new node numbers to old 
node numbers. JNEW (10) gives the new node 
number of the old node number 10. This 
array Is complementary to the array JOLD. 

JOLD (400) 

RENUM 

Array which relates numbers to the new node 
numbers. JOLD (IN) gives the old node 
number of the new node IN. This array is 
complementary to the array JNEW. 

LEM 

CONGF 

Elements used In the G-cal cul at ion * 

LIST (200) 

CL1ST 

Array containing the degrees of freedom 
with a prescribed zero displacement. 

LLIST 

— 

Total length of the LIST array, currently 
LL 1ST - 200. 

LNOD 

— 

Total length of the NOD array, currently 
LNOD - 200 * 8 = 1600. 

LOLD 


Total length of the JNEW and JOLD arrays, 
currently LOLD - 400. 

LX 

- 

Total length of the x array, currently 
LX = 400 * 2 - 800. 

NBAD (200) 

— 

Array which stores the numbers of elements 
that do not satisfy equilibrium conditions. 

NBAND 

INTGR 

Maximum number of rows In the BIGK array. 
Also corresponds to the maximum allowable 
bandwlth, currently NBAND = 75. 



VARIABLE NAME 


COMMON 


DESCRIPTION 


NBOUN 

NBW 

ND (400,10) 

NDA (400) 

ND [S 
NFREE 
NGAUSS 
NCD (4) 

NCK (2) 

NMAT 

NMAX 

NNODE 

NOD (200,8) 


INTGR 


INTGR 


INTGR 

INTGR 

INTGR 

GCAL 


GCAL 


INTGR 


INTGR 


BNOD 


Number of boundary conditions. (Degrees of 
freedom with a prescribed zero displace- 
ment). This value should not exceed LL1ST. 

Current bandwidth of the finite element 
model . 

Array containing the number of nodal ad- 
joining a specific node in each material. 
The first subscript corresponds to a node 
number; the second subscript indicates the 
material number. Array is equivalenced 
with BIGK. 

Array containing the number of nodes 
adjoining a specific node. Array is 
equivalenced with BIGK 

Current degrees of freedom of the finite 
element model, NDIS = NPOIN * NFREE. 

Number of degrees of freedom per node, 
currently NFREE = 3. 

Number of Gaussian integration points in 
each direction. Currently NGAUSS = 3. 

Array containing the node numbers where the 
relative displacements are calculated and 
used in the G-calculation. 

Numbers of nodes at which forces are calcu- 
lated and used in the G-calculation. 

Total number of materials in the problem. 
Currently up to 10 materials can be 
specified. 

Maximum number of columns in the BIGK 
array. Also corresponds to the maximum 
allowable degrees of freedom in the model. 
Currently NMAX = l 110. 

Number of nodes on each element. For the 
parabolic element used in this program, 
NNODE = 8. 

Nodal connectivity of each element. The 
first subscript indicates the element 
number and the second subscript indicates 
the 8 nodes on the element. 
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VARIABLE NAME 

NPLY 

NPOIN 

POUT 

PDER (18,8) 


SHORT 

SMECH 

SP (400,6,10) 

SPA (400,6) 

THETA (10) 
TSTR (8,4) 


C OMMON 

CONFG 

INTGR 

OUT 

CDER 


DESCRIPTION 

Number of pi les used in the laminate. 

Total number of nodes in the finite element 
model . 

Output option Indicator. 

Array containing the parent derivatives, 

9N 9N 

-;r~- at each of the Gaussian points. 

9C 9n 

Array Is stored as, 



j = 1, NGAUSS x NGAUSS. 

OUT Short output option. 

CONFG Mechanical strain £ x prescribed on the 

1 amlnate . 

Average nodal stresses at each node in each 

material due to the combined loading. The 
first subscript Indicates the node number; 
the second subscript indicates one of the 
six stresses o x , o y , o z , c> xy , a , o zx ; 
the third subscript Indicates trie material 
number. Array Is equlvalenced with BIGK. 

Average nodal stresses at each node In the 

model due to combined loading. The first 
subscript indicates the node number; the 
second subscript Indicates one of the six 
stresses o x , o y , o z , o xy , o yz , « zx . Array 

is equivalenced with BIGK. 

MATER Fiber angle for each material. Currently, 

up to 10 materials can be specified. 

CDER 2x2 Gauss lan-Nodal stress transformation 

matrix. 


26 



VARIABLE NAME 

COMMON 

DESCRIPTION 

URL (24,3) 

BLOAD 

Consistent nodal loads of an element. The 
first subscript indicates one of the 
degrees of freedom of the element; the 
second subscript indicates one of the three 
loading conditions of e x = 10~ , AT * 1°, 
or AH = 1%. 

UMO (24,1) 

BLOAD 

Consistent hygroscopic vector at the inte- 
gration point corresponding to AH - 1%. 

UO (24,1) 

BLOAD 

Consistent mechanical vector at the 
integration point cor responding to a 
mechanical strain = 10 

UT (24,0 

BLOAD 

Consistent thermal vector at the 
integration point corresponding to AT = 1°. 

WRIGHT (8,8) 

CAUSS 

Gaussian weights up to 8-point Gaussian 
integrat ion. 

WIDTH 

CON KG 

Width of the laminate. 

X(400, 2) 

BNOD 

Y- and z-coordinate array of all nodes in 
the finite element model. 

XE (8,2) 

— 

Nodal coordinates of the nodes for the 
element being processed. 

XLONG 

OUT 

Long output option. 


i 
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B.3 Common Blocks 


BLOCK NAME ELEMENTS ROUTINE 


AINV 

A(8,8) 

ASEMBL, ASTAR, BMAT 


ASTIF 

BIGK(75, 1 1 10) 

ADJUST, ASEMBL, BOUND, FORCES, 
Q3DC 


BLOAD 

U0(24,l), UT(24, 1) 
UMO(24,l), UEL(24,3) 

BMAT, FORCES, I NT 


BNOD 

X(400,2), NOD(200,8) 

ADJUST, ASEMBL, CKBNDW, CLD, 
GEVAL, PICKF, Q3DC 


CDER 

PDER (18,8), 
TSTR (8,4) 

ASEMBL, FORCES 


CLIST 

. LIST( 200) 

ADJUST, ASEMBL, BOUND, CLD, 
FORCES, Q3DC 


CONFG 

NPLY, IRHS, WIDTH, 
SMECH, DELT, DELM 

ASEMBL, FORCES, GEVAL, PICKD, 
PICKF, STRESS, Q3DC 


DISP 

DIS( 1110,3) 

ASEMBL, BOUND, FORCES, PICKD, Q3DG 


GAUSS 

C0RD(8,8) , WEIGHT(8,8) 

BLOCK DATA, INT, LNTSX, STRESS 


GCAL 

LEM( 2) , NGF(2), NCD 
(4)', GFOR (2,3,3), 
GDIS (2,3,3) 

FORCES, GEVAL, PICKD, PICKF, Q3DC 


GENRL 

GCR(8 , 2) 

ASTAR, BLOCK DATA, FORCES 


INTGR 

NPOIN, NBOUN, NELEM, 
NFREE, NMAX, NBW, 
NBAND, NNODE, NDIS, 
NGAUSS 

ADJUST, ASEMBL, BOUND, CKBNDW, 
CKDEGF, CLD, FORCES, GEVAL, Q3DC 


MATER 

IMAT(200) , EMO( 10,16), 
THETA(IO) 

ASEMBL, FORCES, MODULUS, Q3DC 


MOD 

EE(6>6, 10) , AMOD(6,6), 
ETH(6, 10), ET(6, 1 ), 

IMT, ETM(6, 10) , EM0(6,l) 

ASEMBL, BMAT, FORCES, INTSX, 
MODULUS, STRESS, Q3DC 


OUT 

POUT, SHORT, XLONG 

ASEMBL, BLOCK DATA, FORCES, GEVAL, 
Q3DG 

- 

RENUM 

JOLD(400) , JNEW(400) 

ADJUST, ASEMBL, FORCES, GEVAL, 
PICKD, PICKF, Q3DC 
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APPENDIX C: EXAMPLE PROBLEMS 


In this appendix the input data and detailed outputs for three example 
problems are presented. 

Example Problem 1 

The first example problem Is that of a long rectangular [±35/0^90] g 
laminate with a semi-width b of 0.75 inches. The laminate is edge delami- 
nated at the 0/90 interfaces as Indicated by the upward arrow in the stacking 
sequence. The delamination length a is 0.129 inches. The laminate is sub- 
jected to a uniform axial strain of 0.00254 in/in, a temperature change 

AT of -280°F (cooling from 350°F cure to room temperature of 70°F), and a 
moisture change AH of 0.6%. 

The laminate was modeled with 102 eight-noded elements involving 367 
nodes as shown In fig. 2. The details of the modeling near the delamination 
tip are shown in fig. 5 along with the elements and nodes involved in the 
G-calculat ion . The input data for this model is presented in fig. 6. The 
input data of fig. 6 were obtained using the rectangular mesh generation 
program, DATGEN. This program is described in Appendix E. 

The detailed output from Q3DG for the input in fig. 6 is presented in 
fig. 7. The program calculated the total strain-energy-release rates for M, 

M + T, M f II, and M + T + II loading conditions as 0. 1102, 0.2637, 0.0700, and 
0.1332 lb in/in^, respectively. The CLT analyses of ref. 3 yielded the total 
strain-energy-release rates of 0.1104, 0.2577, 0.0686, and 0.1314 lb in/in^ 
for the same loading cases. The Finite element values, thus, agreed very well 
with CLT values in these cases. 


29 



Example Problem 2 


The second example problem uses the same laminate configuration and load- 
ing as the first problem except that the stacking sequence for this laminate 
is [±45/0/90] o . The input data is the same as shown in fig* 6 except for the 
material specification cards* The material specification cards for this lami- 
nate are presented in fig. 8. The short output option was exercised and the 
output obtained from Q3DG is presented in fig. 9. 

Example Problem 3 

This example is included to demonstrate the input required for the re- 
numbering option. Consider a [0/90] o laminate with edge delaminat Lons between 
the 0° and 90° plies with b = 0-75 in. and a = 0.15 in. The laminate is 
subjected to e x = 0.00254 in/in, AT = -280°F and AH = 0.6%. 

Two identical very coarse finite element idealizations were used. In the 
first model, the nodes were numbered in the thickness direction as shown in 
fig. 10(a). In the second model, the nodes were numbered in the width direc- 
tion as shown in fig. 10(b). Obviously, the width-wise numbering scheme 
yields a larger bandwidth (bandwidth = 59) than the thickness-wise numbering 
scheme (bandwidth - 39). Therefore, the renumbering option was exercised for 
the model shown in fig. 10(b). The input data for the model ;of fig. 10(b) is 
presented in fig. 11. Both models of fig. 10, when executed by the Q3DC 
program, gave identical results for displacements, stresses, and strain- 
energy-release rates. The summary of the strain-energy-release rates for 
these models is shown in fig. 12. 

Typical Results 

Table .1 presents the strain-energy-release rates for 4 edge-del am I nated 
laminates subjected to = 0.001, AT = 280°F and AH = 1.0%. These 
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results ran he used by the us»‘r for verification purposes. These results also 
compare laminates with symmetric edge riel ami nat I ons (about midplane) with 
those with only one edge delamination. 
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Table 1: Strain-energy-release rate results for edge delaminations under lami- 
nate combined loadings: c x = 0.001 ; AT = -280°F; AH = 1.0. Material 

properties used here are: E 1( = 18.2 Msi; E 22 = E 33 = 1.23 Msi; 
v 12 = V 1 3 = v 23 = 0.292 G 12 = G n = G 23 = 0.832 Msi 



Laminate: 

[35/-35/0y|90] s 


M 

M f- T 

M f T h H 

G (Total) 

0.0150 

0.0926 

0.0098 

Gj /G 

94.43% 

43.95% 

88 J\Yl 

C„/G 

5.50% 

56.03% 

11.25% 

G iir /G 

0.07% 

0 . 02 % 

0.30% 

Laminate: [ 35/— 35/0/j^ 90 2 /0/- 

3 5 / 3 5 ] T 


M 

M + T 

M + T f II 

G (Total) 

0.0266 

0.1192 

0.0150 

Gj/G 

75.87% 

36.79% 

97.23% 

G u /G 

24.06% 

63.19% 

2.40% 

G II[/ G 

0 

0 

0 

Laminate: [ 30/-30 2 /30/j90 ] s 


M 

M h T 

M + T + 11 

G (Total) 

0.0273 

0.1292 

0.0110 

Gj/G 

63.58% 

63.58% 

63.58% 

C u /G 

36.90% 

36.90% 

36.90% 

G tII /G 

-0.48% 

-0.48% 

-0.48% 

Laminate: [ 30/-30 2 /30yj90 2 /30/- 

-30 2 /30] T 


M 

M + T 

M + T + 1 ! 

G (Total) 

0.0459 

0.2149 

0.0186 

G z /G 

68.25% 

68.25% 

68.25% 

G ir /C 

32.03% 

32.03% 

32.03% 


Grrr/G 


-0.28% 


-0.28% 


-0.28 % 




APPENDIX D: EXECUTION ERRORS AND DEBUG STRATEGIES 


This appendix presents the conditions under which the program will 
terminate execution and discusses potential debug strategies. 

The program will terminate execution and print a self-explanatory 
diagnostic message under the Following conditions: 

1. A node or nodes have unspecified coordinates. 

2. An element has an unspecified material group number. 

3. The degress of freedom In the model is greater than NMAX, the column 
dimension of BIGK. 

4. The bandwidth of the model Is greater than NBAND, the row dimension 
of BIGK. 

5. Negative diagonal terms are computed in the element stiffness matrix 
of any element. 

6. The band-solver aborts the program under the following conditions: 

a. The matrix is not positive definite. 

b. The bandwidth is greater than the number of equations. 

The first two error conditions can be corrected by checking and modifying 
the Input. The third and fourth error conditions are caused by the size of 
the model. To correct these errors, two alternatives are available. The 
first Is to increase the dimensions of BIGK and the other associated variables 
of the program, as explained earlier. The second approach is to reduce the 
size of the finite element model. Error condition 5 is usually caused by an 
improper description of one or more elements. Each element's nodal connect- 
ivity must be described In the counterc lockwi se sense sta rti ng at any co rner 
node^ If this is not done correctly error condition 5 may result. Another 
possible source of this error may be erroneous nodal coordinates. The band- 
solver's error condition a Is due to failure to specify sufficient restraints 
to the model to prevent rigid body motion. Therefore, the boundary conditions 



specified on the model should be checked and corrected. Error condition b of 
the band-solver results whenever the bandwidth is greater than the total 
degrees of freedom in the model. This condition usually occurs because of 
erroneous nodal numbering. This condition, however, will never occur because 
the program would abort execution before it calls the band-solver. 

Plotting the finite element model is recommended since a plot will 
quickly reveal any input errors in the nodal coordinates and element connec- 
tivities, which constitute the bulk of the input. Any other errors which are 
encountered can usually be debugged using the host computer debug and dump 
utilities. 

As an additional check, the total strain-energy-release rates computed by 
Q30G should be compared with those obtained by CLT theory analyses [1-3] which 
can be performed on microcomputers. The Q3DG results should agree very well 
with the CLT results even for coarse models. Therefore, if there is a large 
discrepancy between the Q3DG results and CLT results, both analyses should be 
examined and the causes for these discrepancies should be determined. 

! 
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DATGEN PROGRAM 


APPENDIX K: 

This appendix describes a rectangular mesh generation program, DATGEN, 
which creates a data file that can be used as an Lnput for the Q3DG program. 
First, the capabilities of this program are explained. Next, a few examples 
of data generation and types of errors encountered in the use of DATGEN are 
presented. Finally, the input data for the program is given. 

E.l Capabilities of DATGEN 

DATGEN is a program which can be executed interactively on any mainframe 
computer. The program requires 41,000 octal words of core memory with the 
current dimensions. With the Input of y- and z-coordinates on the 
y 3 constant and z - constant lines and the location of the delamination tip, 
the program generates a rectangular mesh, computes the nodal coordinates of 
the midside nodes and the nodal connectivities. The user then enters the rest 
of the data requested by the program. The program formats the input and at 
the end of the interactive session, a complete data file is ready for use as 
an input file for Q3DG. 

Material properties and thermal expansion coefficients of the T300/5208 
graphite/epoxy unidirectional composite lamina are built into the DATGEN pro- 
gram’s BLOCK DATA. The user can change these properties very easily using a 
text editor. ' 

The input variables of DATGEN are listed in section E.4. 

E.2 Examples of Interactive Sessions 

Consider the [±35/0^90] s laminate of example 1 in Appendix C. A finite 
element model for this laminate was shown in fig. 2. The y- and z-coordinates 
of y - constant and z = constant lines used in this model are presented in 
fig. 13 and were stored in file YZ. There are 18 y-coordinates and 7 
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z-coordinates listed In this file which correspond to the y- and z-coord inatos 
of the vertical and horizontal lines, respectively. Therefore, the values of 
NYP and NZP in the program are 18 and 7, respectively. The delamination tip 
is located at (0.6210, 0.0054). The y- and z-coordinates of the delamination 
tip correspond to the llth y-coordinate and the 3rd z-coordinate in the file 
YZ. Therefore, the values NYC and NZC in the program are 11 and 3, 
respectively. 

Fig. 14 presents the complete interactive session for this model. All 
input is echoed onto the terminal screen. The output from DATGEN is written 
on the file TRIAL. A complete listing of TRIAL is presented in fig. 6. 

E.3 Recoverable and Irrecoverable Errors 

During an interactive session it is possible that the data may be input 
incorrectly. Some errors in the data input can lead to recoverable errors 
while certain others can lead to irrecoverable errors. This section presents 
examples of each kind. 

Fig. 15 presents a partial interactive session for the [±35/0 / 9 0 1 _ 

t s 

laminate discussed in fig. 14, where errors were made during the input of 
material and thickness specifications. The bottom line of the ply did not 
correspond to any of the given z-coordinates. As illustrated by fig. 15, the 
errors are recoverable. The program warns of improper Input and recycles back 
to the point where the error was made and allows the user to enter the correct 
data. 

Fig. 16 illustrates an interactive session that terminated due to an 
irrecoverable error. If the sum of the thicknesses of all plies does not 
equal to the maximum z-coordinate input by the user, DATGEN terminates the 
execution with an error as shown in fig. 16. 
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E.4 Renumbering Option 


I*Mg. 17 presents an Interactive session in which a finite element model, 
as shown in fig. 10(b), is created for a [0^90] g laminate with a delamination 
between the 0° and 90° plies. In this session the renumbering option was 
exercised. The renumbering scheme was stored in the file REN and the y- and 
z-coordinate were stored in the file YZC. These two files are shown in 
fig. 17. 

When the renumbering option is exercised, DATGEN checks the renumbering 
scheme to see if any number is repeated, exceeds the maximum number of nodes 
in the model, or does not occur at all. If any of the three conditions 
occurs, DATGEN gives a summary of the error condition and no renumbering 
scheme Is entered on the output data file. An example of DATGEN* s warning 
message for this error Is presented In fig. 18. In this example, a node with 
number 39 was repeated twice and a node with number 30 was not given in the 
renumbering scheme. 
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E.5 Description of Input Data for DATGEN 


CARD SET 

COLUMNS 

FORMAT 

VARIABLE 

DESCRIPTION 

1. 

1-10 

A10 

CORCON 

Name of flip on which output is 
written. This Is tin 1 f 1 1 < ‘ that 
will be used as Input for Q3DC. 

2. 

1-80 

8A10 

TITLE 

Title of the problem. 

3. 


* 

NYP, NZP 

Number of the corner nodes in y- 
and z-dl rec t ions , respectively, of 
the rectangul ar mesh. 

4. 


* 

NYC, NZC 

Number of the y- and z-coordinates, 
respectively, of the delamination 
tip. 

5. 

1-10 

A10 

CFILE 

File name on which the y- and 
z-coordinates values are stored. 

Note: 2Z. 

and z-coordinates should 

be greater than or equal to zero . 

6. 


* 

NMAT 

Number of materials in the model 
(each ply with a different fiber 
angle is classified as a different 
material ) . 

7. 


* 

THETA 

Fiber angle for each material. 

This line 

occurs once 

for each material and 

hence occurs NMAT times. 

8. 


* 

NPLY 

Number of piles in the laminate. 

9. 


* 

MATC 

Material number of each ply. 

10. 


* 

PLH 

Ply thickness. 

The lines 

9 and 10 occur NPLY times* The p 

rograrn assumes that the plies are 

numbered from the 

top of the laminate. 


ii. 


* 

IRHS 

Number of loading conditions. Note 


that the number of loading 
conditions IRHS < 3. 


Note: * denotes a free format* 
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CARD SET 

COLUMNS 

FORMAT 

VARIABLE 

DESCRIPTION 

12. 



* 

SMECH 

Mechanical strain e x . 

13. 



* 

DELT 

Temperature change, AT, 

14. 



* 

DELH 

Moisture (hygroscopic) absorption, 
AH. 

If IRHS = 1, 
needed. 

Card 

Also 

set 

note 

13 and 14 are not needed, 
that Input of e x = 0 is 

If IRHS = 2 Card set 14 is not 
not permitted by Q3DG. 

15. 

* 


* 

IRENUM 

Renumbering option. IRENUM = 0, no 
renumbering is attempted. 

IRENUM = 1, renumbering scheme 
follows. 

16. 

1-10 


A10 

RF LIE 

Name of file which contains the 
renumbering scheme. If a blank is 
input, renumbering scheme should be 
input manually. 

17. 

* 


* 

JNEW (1) 

• • • 

JNEW (NPOIN) 

Renumbering scheme manual input. 
Separate numbers by commas, use 
as many lines as required. 


Note: * denotes a free format 
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Spe cial Cases: 


The above input will generate the rectangular mesh and data file for 
delarainat ions at any location other than the midplane. When the delaminatLon 
is at the midplane, input NYC = NYP and NZC = 1. The resulting data is 
correct except the G-calculat ion cards, card sets l l, 12, and 13, and tin? 
boundary conditions. The user can use the text editor and input the correct 
integers in card sets 11, 12, and 13 and delete the unnecessary boundary 
conditions. Also, be sure to use the negative integers on card set 13 as 
explained in the special cases part of the INPUT DATA section. 

If the user wants to perform only a Q3D edge-stress analysis, then as 
♦ 

above, use NYC *= NYP and NZC = 1. Furthermore, change card sets 11, 12, 
and 13 to zero values as explained in the special cases part of the INPUT DATA 
section. The boundary conditions generated by the DATGFN are correct for this 
case. 


APPENDIX F: COSMIC DOCUMENTATION 


This appendix describes the program documentation requirements required 
by COSMIC. Items such as method of solution, user instructions, accuracy of 
results, sample input and output, and flow chart were already discussed in 
detail in the text and Appendixes A-E. Therefore, this appendix limits itself 
to the remaining items described in COSMIC software submittal guidelines. 

1. Computer Configuration : The Q3DG program was written, compiled and 

successfully executed on CDC Computers such as on CYBER 173, 175, 855 
using NOS 1.4 and 2.1 operating systems. 

2. M emory R eq uired : With the current dimensions Ln the program, about 

343,000 octal words of memory were needed for execution. As pointed 
out earlier, with the current dimensions, finite elements models of up 
to 370 nodes with a bandwidth of 75 can be successfully run with the 
program. 

3. Source Language : The Q3DG program was written in FORTRAN V. The 

compilation showed that, except for one statement in BLOCK DATA 
subprogram (where the output options are set), the program is machine 
independent FORTRAN. This statement DATA SHORT, XLONG/ 5HSH0RT, 
5HXL0NG/ can be changed, for example, to DATA SHORT, XLONG/ 1 111. 11, 
9999.99/. However, this was not done because SHORT and LONG clearly 
are meaningful to describe the output option. 

4. Program Timing: Execution times varied widely and are dependent on 

the size of the finite element models. Small models with less than 80 
nodes were executed within about 20 CP seconds interactively. The 
largest model with 367 nodes presented in Appendix C took about 184 CP 
seconds on CYBER 855 Computer. 
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5. T ape Documentat io n: The tape containing the Q3DG program Is a 

9 track, 1600 BPI, CDC NOS Internal unlabeled tape containing eight 
files. The first file contains the source code of the Q3DG program. 
The second file contains the source code for the rectangular mesh and 
data generator DATGEN. The third file contains input data for the 
[±35/0/90] laminate, shown in fig. 6, discussed in Appendix C. The 

T S 

fourth file contains the output of Q3DG for the input shown in 
fig. 6. This output file is presented in fig. 7. 

The fifth file of the tape is the input data for the [±45/0/|90] s 
laminate, discussed in Appendix C. The sixth file contains the output 
of Q3DG for the input data in the fifth file. This output file is 
presented in fig. 9. 

The seventh file of the tape contains the y- and z-coordinate 
file shown in fig. 13. This file is used in conjunction with the 
DATGEN program and the resulting interactive session is the last 
(eighth) file of the tape. 

The NOS procedural file which created this tape is shown in fig. 19. A 
microfiche which contains the compiler listings of Q3DG and DATGEN, the input 
and output files, and interactive sessions of the DATGEN program is enclosed 
with this manual. 

Finally, a dayfile of a run on the NOS 2.1 operating system is presented 
in fig. 20. The output is written and saved on a file SOUT. The dayfile is 
written and saved on a file DAY. The dayfile is also included in this figure. 
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APPENDIX G: PC FLOPPY DISKS 


Because of the popularity of the IBM PC and its compatibles, the Q3DG 
program was downloaded on to floppy disks. The floppies titled Q3D-1 and 
Q3D-2 are enclosed with this manual. These disks contain the files listed 
under Tape Documentation of Appendix F (see page 42). The directories for 
these files are reproduced below. 

Any one of the PC communications packages, such as PC-PLOT , PC-TALK , 
etc., can be used to upload all the files from the Q3D floppy disks on to the 
mainframe computers easily. The procedures to do this can be found in the 
documentation of the communication packages. 


A: \>DIR/P 


Volume in drive A is Q3D-1 
Directory of A:\ 


Q3DG 

DATGEN 

TFICH1 


82688 
17664 
240768 
3 File (s) 


9-05-86 2: lip 

9-09-86 10:03a 

9-08-86 3 : 16p 

19456 bytes free 


A: \>B: 


B: \>DIR/P 


Volume in drive B is Q3D-2 
Directory of B:\ 


TFICH2 

TFICH3 

TFICH4 

TFICH5 


245248 

32128 

61312 

8064 

4 File(s) 


9-08-86 4 : OOp 

9-08-86 4 : 13p 

9-09-86 9:34a 

9-09-86 9:55a 

14336 bytes free 
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APPENDIX H: LIST OF SYMBOLS 



G » G I » G i I » G m 


h 

AT 


t 

U,V,W 


u,v,w 


x,y ,z 


a. 

1 

e 

o 

A 

0 


{ 0 } 

{e} 


delarainat Ion length 

ha If -width of the laminate 

stress-strain matrix 

Young's modulus in the i-direction 

force at node i in the x~, y-, and z- directions* 
respect Lvely 

total, mode I, mode II, and mode Hi stra in-energy release? 
rates , respectively 

shear modulus 

ply thickness 

temperature change 

half-thickness of the laminate 

displacement functions 

displacements in x-, y~, and z-direct ions , respectively 
Cartesian coordinates 

coefficient of thermal expansion in the i-direction 
uniform axial strain imposed on the laminate 
length of the elements nearest to the del ami nation tip 
ply fiber angle measured from the x-axis toward the y-axis 
Poisson’s ratio 

Cartsian stresses, {o o o a a o } 

x y z xy yz zx 

Cartsian strains, {e e e e e e } 

x y z xy yz zx 


Subscript 


i,j i,j = 1>2,3 

1,2,3 longitudinal, transverse, and thickness directions of a zero 

degree ply 
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Fig. 2 A Typical Finite Element Idealization 




Fig. 4 Flow Chart for Q3DG Program 
(Read the flow chart from top to bottom and left to right) 
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Fig. 5 Detail Near the Delamination Tip for the Model Shown 
in Fig. 2 and for a Midplane Delamination 
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Fig. 7 (Continued) 
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Fig. 7 (Continued) 
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Fig. 7 (Continued) 
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Fig. 7 (Continued) 


SOME MORE PROBLEM SPEC I F IC AT IONS 
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Fig. 7 (Continued) 
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Fig. 7 (Continued) 
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Fig. 7 (continued) 
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(Concluded) 








1950000E+08 . HP??.- .'C+t>7 . 14B0000E+07 . 3000000E+00 


U) OJ G «d -* OJ IS tO *4 OJ G uo -4* OJ 
SQSQQSGSiSQSfJSQS 
4 114 4 1 14 4 1 I 4 4 I I 

LU Id Id Id L! Id UJ UJ UJ Id UJ U Id til Id 

GGGGGGGGG G G G G G ® 
| SSi < S , S^O , S , S , S6i ,, S^t<i G G 
GGGGGGGGGGGGGG® 
S<SG< 5 GG'&^GGG'SSiSQ 
G(7iidGGC*>OGG C^<J3 Q <S (J'O) 
G ^ ID G G <4 f) Q 5» 4- in g g 4 - in 

cd ^ m io <o -« in n co ^ in n co w in 


vd id OJ 

N 

id id 

•M 

r* id 

^ OJ N 

id 

id OJ 

G G G 

G 

G G 

G 

G G 

G G G 

G 

G G 

+ l 1 

4 

4 l 

l 

4 4 

l 1 4 

4 

1 t 

Id UJ UJ 

Id 

U id 

Id 

Id Id 

IU Id tu 

l;J 

111 LU 

G G G 

G 

G G 

G 

G G 

G G G 

G 

G G 

G G G 

G 

G G 

G 

G G 

G G G 

G 

G G 

G G G 

G 

G G 

G 

G G 

G G G 

G 

G Q 

G G G 

G 

S G 

G 

G G 

G G G 

G 

G G 

G G id 

<0 

S G 

id 

CO G 

G uJ 03 

G 

G id 

G M in 

* 

G M 

in 

^ G 

fO tn 

G 

ro Ul 

<o oj in 

• • • 

• 

CO OJ 

• • 

m 

9 

H CO 

• • 

co m *-< 

• • • 

CO 

• 

cu ui 

• • 


I I 


G i* J G N 

G OJ 

G N G 

© G N 

G 

OJ G 

G G G 

G G G 

G G G 

G G G 

G 

O G 

+ 

9- 

4 

4 

* 

4 

4 + 4 

4 

4 

4 

4 

4 4 

’ J 

U 

IJ 

i,* 

; j 

i .; 

V ;.» !,; 

U 

UJ 

l:l 

• il 

UJ UJ 

r . 



f w 

* , 

/ # 

* , G * , 

* 

G 

G 

G 

G G 

*, 

• , 




* . 

' • . ' # 

‘ , 

* J 

G 

’ * 

G G 

G 

G 

G 


■* V 

* < 

G G G 

G 

G 

G 

G 

G G 

G 

G 

i» j 



t ; 

* . f i n 

G 

G 

G 

G 

G © 

G 

G 

G 

C 

r t 

G 

* * C G 

G 

G 

CO 

G 

G O 


L’l 

G 


f % 

L ’ 

G J G 

G 

G 


G 

G G 

'•'i 

«? 

G 


r) 

d 

G N 

G 

G 

**• 

ri 

O 


I 


© vd 

(0 

G 

id 

* 

CO 

G 

u3 

-t co 

G 

uD >t 

G G ® 

G 

G 

& 

G 

G 

G 

G 

G G 

G 

G G 

4 4 1 

4 

4 

4 

1 

4 

4 

4 

1 4 

4 

4 1 

UJ UJ UJ 

UJ 

UJ 

UJ 

UJ 

UJ 

UJ 

UJ 

Id UJ 

tu 

U LU 

G G & 

IS 

G 

G 

G 

G 

G 

G 

G G 

G 

G G 

& G G 

G 

G 

G 

G 

G 

G 

G 

G G 

G 

G G *-« 

G G G 

G 

G 

O 

G 

G 

G 

G 

G G 

G 

G G 

G G & 

G 

G 

O 

S 

G 

G 

G 

G G 

G 

G G 

G & O' 

in 

G 

G 

CO 

U1 

G 

G 

C 7> in 

G 

G td 

G G *4* 

c* 

G 

G 

>4 

CU Q 

G 

•>4 (d 

G 

G 

W CO 

t-4 

10 cO 


^ co 

<0 

*-« 

n 

CO ^ OJ 


G 


* 


*0 

<0 


CO 

cn 

•rl 

Ll_ 


I 


77 


Partial listing of input data for 
Example Problem 2 . 
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Fig. 9 (continued) 
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Fig. 10 Very Coarse Models for a [0/90] Laminate 
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Input data for the model shown in fig. 10 (b), . 
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Fig. 11 (Concluded) 
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Fig. 12 Strain energy release rates obtained with the 
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Fig. 12 (Concluded) 
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Fig. 13 Listing of the coordinates file, YZ. 
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Fig. 15 (Concluded) 
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Fig. 16 (Concluded) 
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Fig. 17 (Concluded) 
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(Concluded) 



RAJU 


Q3DGT , T200 . B 1205 

USER, User number, Password. 

CHARGE, Charge number, LRC. 

GET , Q3DG, DATGEN , D3 5 , D4 509 , YZ , TINT , SOUT, SOUT2 . 
REWIND, *. 

LABEL ( TNAME , NT , PO=W , LB=KU , F=S , VSN=Q3 DG ) 

TCOPY ( Q3 DG , TNAME , F=E , CC=8 0 ) 

WRITEF, TNAME. 

TCOPY ( DATGEN , TNAME , F=E , CC=8 0 ) 

WRITEF, TNAME. 

TCOPY ( D3 5 , TNAME , F=E , CC=8 0 ) 

WRITEF, TNAME. 

TCOPY ( SOUT , TNAME , F=E , CC=1 3 6 ) 

WRITEF, TNAME. 

TCOPY ( D4 5 0 9 , TNAME , F=E , CC=8 0 ) 

WRITEF, TNAME. 

TCOPY (SOUT2 , TNAME , F=E , CC=13 6 ) 

WRITEF , TNAME . 

TCOPY ( YZ , TNAME , F=E , CC=8 0 ) 

WRITEF, TNAME. 

TCOPY ( TINT , TNAME , F=E , CC=1 3 6 ) 

WRITEF, TNAME. 

DAY FILE, DT APE . 

REPLACE , DTAPE . 

EXIT. 

DAYFILE, DTAPE. 

REPLACE, DTAPE. 


FIG. 19 NOS Procedural file used to create the C0SH1C t’ape. 
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FULL DAYFILE. 85/08/06. 20.12.09. 

JOB NAME A3MIHUM ON C-MACHINE. 

09.07.04. LQ6 , T300 . B 1205 

09.07.04. RAJU 

09 . 07 . 04 . USER , 654237E , . 

09.07.04. CHARGE , 1 01 955 , LRC . 

09.07.05. GET , 03DG , TAPE5=D35 . 

09.08.49. FTN5 , I=Q3DG , L=OUT . 

09.12.39. 62700 CM STORAGE USED. 

09.12.39. 33.685 CP SECONDS COMPILATION TIME. 

09. 12.40. REWIND, + . 

09.12.40. 6 FILE(S) PROCESSED. 

09. 12.40. DELIVER. B 1205 RAJU. 

09.12.40. ROUTE , OUT , DC=LP . 

09.12.40. ROUTE COMPLETE. 

09.12.42. LDSET ( PRESETA=NGINF , MAP=SBEX ) 

09.12.42. LG0(TAPE5,S0UT) 

20.12.06. STOP 

20.12.06. 343600 MAXIMUM EXECUTION FL. 

20.12.06. 185.084 CP SECONDS EXECUTION TIME. 

20. 12.06. REWIND,*. 

20.12.06. 7 FILE(S) PROCESSED. 

20.12.06. REPLACE , SOUT . 

20. 12. 09. DELIVER. B 1205 RAJU. 

20.12.09. ROUTE , SOUT , DC=LP . 

20.12.09. ROUTE COMPLETE. 

20. 12. 09. DAYFILE, DAY. 


FIG. 20 A sample dayfile of a NOS run of Q3DG. 
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